######################################################################
#Replication file for "The Impact of China's AIIB on the World Bank"
#Authors: Jing Qian, James Vreeland, Jianzhi Zhao
#Codes to replicate Figure A.4
######################################################################
rm(list = ls())

#To install the bpCausal package (ver 0.0.1)
# install.packages('devtools', repos = 'http://cran.us.r-project.org') # if not already installed
# devtools::install_github('liulch/bpCausal')

#Load bpCausal package
library(bpCausal)

#Other packages used
library(tidyverse)
library(RColorBrewer)

#Setwd
setwd("C:/Users/qianj/Dropbox (Princeton)/AIIB_WB_Replication")

#Load custom functions
#Note: Given the number of tests performed in the appendix, 
#several custom functions are created to implement estimation methods in a more efficient way.
#These functions merely create a wrap-up of existing functions, like bpCausal and gsynth,
#in order to avoid typing parameters repeatedly.

source("code/custom_functions.R")

###################################
#Estimation for Figure A.4
###################################
#----------------
#(1) Load data
#----------------
load("data/data_main.RData")

#--------------
#(2) Setup
#--------------
D.use = "aiib_founder_2016"
dv.list = c("hard_project_count_all")
covs.all = c("gdppc_log_lag",
             "population_log_lag",
             "election_either_lag",
             "fdi_gdp_lag",
             "debt_gni_lag",
             "oda_gni_lag",
             "polity2_lag",
             "unsc_lag",
             "IdealPointDistance_lag")

df.use = df.main

#Container
result.a4 = list()

#--------------
#(3) Estimation
#~22 minutes
#--------------
now = Sys.time()
for (i in 1:length(covs.all)){
  covs.use = covs.all[-i]
  X.use = Z.use = A.use = covs.use
  result.a4[[i]] = bpcausal.group(dv.list = dv.list,
                                  df.use = df.use,
                                  D.use = D.use,
                                  X.use = X.use,
                                  Z.use = Z.use,
                                  A.use = A.use)
  print(i)
}
print(difftime(Sys.time(), now))

#Names
names(result.a4) = covs.all

#--------------
#(4) Number of observations
#--------------
obs.a4 = c()
for (i in 1:length(covs.all)){
  obs.a4[i] = df.main %>%
    .[, c(D.use, dv.list, covs.all[-i], "iso2c", "year")] %>%
    .[complete.cases(.),] %>%
    nrow()
  cat("Exclude", covs.all[i], "N =", obs.a4[i], "\n")
}

###################################
#Produce Figure A.4
###################################
#-------------
#(1) Prepare results
#-------------
result.nocov.use = lapply(result.a4,
                          function(x) effSummary(x$hard_project_count_all))

#Data
est.use = result.nocov.use %>%
  lapply(., function(x) x$est.avg)

coef.use = est.use %>%
  lapply(., function(x) x[1]) %>%
  unlist()

ci.l.use = est.use %>%
  lapply(., function(x) x[2]) %>%
  unlist()

ci.u.use = est.use %>%
  lapply(., function(x) x[3]) %>%
  unlist()

#-------------
#(2) Parameters
#-------------
#Labels
label.use = c("GDP pc",
              "Population",
              "Election",
              "FDI",
              "Debt",
              "ODA",
              "Polity",
              "UNSC",
              "UNGA") %>%
  paste0(., 
         "\n(N=", 
         obs.a4,
         ")")

#Parameters
n = 9
xlim = c(-2, 2)
cex = 2
xtext = 0.6
colors = brewer.pal(n = 3,
                    name = "Dark2")

#------------------
#(3) Plot
#------------------

#------------------
#Setup
png(filename = "figure/Figure_A4.png",
    width = 800,
    height = 1200)

par(mar = c(5, 16, 5, 1))

#----------------
#Empty plot
plot(1, 
     type = "n",
     xlab = "Estimated Average Treatment Effect on the Treated",
     ylab = "",
     yaxt = "n",
     xaxt = "n",
     xlim = xlim,
     ylim = c(n + 0.5, 0.5),
     cex.lab = 2)
title(main = "World Bank Infrastructure Projects",
      line = 1,
      cex.main = cex)
box()

#Points
points(x = coef.use,
       y = 1:n,
       col = colors[2],
       pch = 19,
       cex = cex)
#Arrows
arrows(x0 = ci.l.use,
       x1 = ci.u.use,
       y0 = 1:n,
       y1 = 1:n,
       length = 0.05,
       angle = 90,
       code = 3,
       col = colors[2],
       cex = cex,
       lwd = cex)

#Other elements
#abline
abline(v = 0,
       lty = "dashed",
       lwd = cex,
       col = "gray")
#Axis
axis(side = 1,
     at = seq(min(xlim),
              max(xlim),
              1),
     cex.axis = cex)
#Removed countries
axis(side = 2,
     at = 1:n,
     labels = label.use,
     xpd = NA,
     cex.axis = 2,
     las = 2,
     tck = 0.000001)
#-----------
#Close
dev.off()